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Abstract 

We propose a method to find the QCD critical point at finite density calculating the canon- 
ical partition function Xc{T,N) by Monte-Carlo simulations of lattice QCD, and analyze data 
obtained by a simulation with two-flavor p4-improved staggered quarks with pion mass m„ x 
770MeV. It is found that the shape of an effective potential changes gradually as the temperature 
decreases and a first order phase transition appears in the low temperature and high density re- 
gion. This result strongly suggests the existence of the critical point in the {T,fig) phase diagram. 



1. First order phase transition and canonical partition function 

The critical point terminating a first order phase transition line in the phase diagram of QCD 
at high temperature and density is one of the most characteristic features that may be discovered 
in heavy-ion collision experiments. To understand the phase structure, first principle calculations 
of QCD by numerical simulations are very important. One of the interesting approaches to find 
a first order phase transition is to construct the canonical partition function Xc(T,N) by fixing 
the total quark number (A^) or quark number density (p). From the canonical partition function, 
one can estimate the quark number giving the largest contribution to the grand partition function 
Z,GciT,Hq)- Because two different states coexist at a first order transition point, two different 
quark numbers give equally large contributions simultaneously if the transition is of first order. 

The canonical partition function is defined by a fugacity expansion of XGc(T,Hq), 



cc{T,l^c) = r !DU[AttM{^iq|T)f e-'^ = ^ Zc{T,N)e'"'^l\ (1) 



where detM, S g and A^t are the quark determinant, the gauge action and the number of flavor, 
respectively. The term Zc(T,N)e'^^'i^^ can be regarded as the probability distribution of the 
quark number A^. Moreover, it is worth introducing an effective potential Veff as a function of A^, 

VeffCM T,^lq) = - lnZc(T,N) -N^^ - N^, Zcc(T,fig) = ^ e"^-, (2) 

N 

where / is the Helmholtz free energy. In a first order phase transition region, this effective 
potential has minima at more than one value of A^. At the minima, the derivative of Veff vanishes: 

dN ^ ' dN T ^ ' 
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Figure 1: Phase structure in the (T,p) plane and the behavior of /ig/T as a function of p. 



Hence, in the first order transition region, we expect d(lnXc)/dN(T,N) = -j-tq/T takes the same 
value at different A^. Here, jj.*{T,N) is the chemical potential which gives the effective potential 
at a minimum point (T, N) and becomes fig in the thermodynamic limit. 

The phase structure in the {T,p) plane is sketched in the left panel of Fig. [T] The thick line is 
the phase transition line. We expect that the transition is crossover at low density and becomes 
of first order at high density. In the first order transition region, the two coexisting states are 
mixed. The region between the two thick lines is the mixed state. The expected behavior of fi* 
along the lines A and B are shown in the right figure. When the temperature is higher than the 
temperature at the critical point Tp^ (line A), jj.* increases monotonically as the density increases. 
However, for the case below T^p (line B), this line crosses the mixed state. Corresponding to the 
double-well potential in a finite volume, is expected to be an S-shaped function. In the infinite 
volume limit, fi* does not increase in the region between pi and p2, since the surface energy 
between the two states can be ignored. 

The Glasgow method ^ has been a well-known method to compute the canonical partition 
function. A few years ago, the above mentioned behavior at a first order phase transition was 
observed in 4 flavor QCD with staggered fermions calculating the quark determinant by the 
Glasgow algorithm on a small lattice Also, simulations with a fixed quark number, i.e. in the 
canonical ensemble, have been tried fj]. However, the calculations by the Glasgow algorithm 
and the simulations for a canonical ensemble requke large computational cost and are difficult 
except on a small lattice with present day computer resources. In this report, we propose a 
method based on a saddle point approximation [4]. By this approximation, the computational 
cost is drastically reduced and the first order like behavior was observed for 2 flavor QCD. 

2. Inverse Laplace transformation with a saddle point approximation 

From Eq. ([TJ, the canonical partition function can be obtained by an inverse Laplace trans- 
formation, 

Zc(T,N) = |- e-^^»'^-'^''^^Zoc(T,^io + iHi)d(^) , (4) 

where /uq is an appropriate real constant and /i,- is a real variable. We have used the fact that 
Zcc(T,Hq + 2mT/3) = Xgc(T , /J-cj) for any complex fj.^. The grand partition function can be 
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evaluated by the following expectation value at - 0. 



ZGc{T,Hq) 

Zgc(T,0) "Z, 



GC J 



DU 



detM(//,/r) \^ 
detM(O) / 



detM(0)™'e 



detM(nq/T) 
detM(O) 



■(5) 



We apply a saddle point approximation to evaluate Eq. (|4|i, which reduces the computational 
cost. If one selects the yUo at a saddle point in Eq. dU, the necessary information is limited that 
around the saddle point when the volume is sufficiently large. Moreover, if we restrict ourselves 
to study the low density region, the value of detM(fj.q/T) near the saddle point can be estimated 
by a Taylor expansion around fig = 0. The calculation by the Taylor expansion is much cheaper 
than the exact calculation and the study using a large lattice is possible. Also, the truncation error 
can be systematically controlled by increasing the number of the expansion coefficients. 

We assume the existence of a saddle point zq in the complex fig/T = z plane for each con- 
figuration, which satisfies D'(zo) - p - 0. Here, the quark number density in a lattice unit and a 
physical unit are p = N/N^. and p/T^ - pNf, respectively, (detM(z)/ detM(0))'^f = exp[A^^D(z)] 
and D'(z) = dD{z)/dz- We then perform a Taylor expansion around the saddle point and obtain 
the canonical partition function. 



Zc(T,pV) 



In 



ZGc(r,0) 



/ r exp 



y\D(z^)--pZi)--D"{zo)x' 



dx] 



V2^ 



ZcciT, 0) exp [V (D(zo) - pzo)] e 



-ial2 



1 



V|Z)"(zo)l 



(6) 



where D"{z) - d^D{z)/dz^, V = and D"(z) = \D"{z)\e'". Higher order terms in the expansion 
of D{z) are negligible when the volume V is sufficiently large. 

Within the framework of the saddle point approximation, the derivative of In Zc with respect 
to or p can be evaluated by 



_ 1 d\nZc(T,pV) 
T ~ V dp 



zo exp [V(D(zq) - pzo)] e-'"!^ 



(r,/i,=0) 



exp [V (Dizo) - pzo)] e-'-l^ ^^^5 



(7) 



(r,A',=0) 



This formula is similar to that of the reweighting method for finite ^q. The operator in the 
denominator corresponds to a reweighting factor, and /i*/r is an expectation value of the saddle 
point calculated with this modification factor 



3. Numerical results and conclusions 

We compute the derivative of In Zc using the data obtained in with the 2 flavor p4- 
improved staggered quark action, m„ x 770MeV. Because the operators in Eq. O are com- 
plex, this calculation suffers from the sign problem. To eliminate the sign problem, the ap- 
proximation proposed in [6] is used: If one assumes that the distribution of the complex phase 
is well-approximated by a Gaussian function, the complex phase factor e"' can be replaced by 
exp[-(0^)/2]. We estimate the quark determinant by the Taylor expansion up to 0{fiq). Because 
the calculation of Eq. (|7]) is similar to the calculation by the reweighting method, the configu- 
rations which give important contribution are changed by the modification (reweighting) factor. 
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Figure 2: Chemical potential vs. quark number density for Nf = 2 with a saddle point approximation. 

To avoid this problem, we use the multi-/? reweighting method. By this method, the important 
configurations are automatically selected among all configurations generated at many simulation 
points of iT,fig - 0). The details are given in [4]. 

The result of fJ.*g/T is shown in Fig.|2]as a function of p/T^ for each temperature T/Tc- Tc is 
the pseudo-critical temperature at = 0. The dot-dashed line is the value of the free quark-gluon 
gas in the continuum theory, p/T^ = NiH^qjT) + (l/7r^)(jUq/T)^]. From this figure, we find that 
a qualitative feature of fJ.*/T changes around T/Tc ~ 0.8, i.e. fi*/T increases monotonically as p 
increases above 0.8, whereas it shows an S-shape below 0.8. This S-shape is a signature of a first 
order phase transition. With some approximations, the critical value of f^*/T is is estimated to be 
about 2.4, which is roughly consistent with the critical point estimated in 14^ by calculating the 
effective potential of the plaquette using the same configurations, (T/Tc, //^/T) ^ (0.76,2.5). The 
difference between these two results may include a systematic error. Our result strongly suggests 
the existence of the critical point terminating the first order phase transition line in the iT,fiq) 
phase diagram of QCD. Further studies are necessary to predict the critical point quantitatively. 
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